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The least squares fit to a straight line, when both variables are affected by all equal uncorrelated 
errors, leads to very simple results for both the estimated parameters and their standard errors, of 
widespread applicability. In this paper several formulas are derived, presenting a full set of results 
about the estimated parameters and their standard errors. All the results have been validated with 
extensive Monte Carlo simulations. The emphasis of the paper is on the calculation and properties 
of the best-fit parameters and their standard errors. 



2 



I. INTRODUCTION 



Least squares fit to a straight lin e (LSFSL) is the subject of an extensive hterature, not only in the field of physics. 
See, for instance, the reviews and the references therein to the older literature. 

A good understanding of the LSFSL is important, not only for students but also for researchers in physical sci- 
ences, because least squares fits, nowadays, are typically done by black-box computer programs. On the other hand, 
the availability of fast computer programs allows one to carry on long and intensive calculations and Monte Carlo 
simulations in a short time, often avoiding the use of approximations. 

In physics, problems of LSFSL are often met such that both variables are affected by significant measurements errors, 
equal for all measured data points separately for both variables, and with the errors in two variables uncorrelated; 
this is the so-called Standard Weighting ModeP^ (SWM). In fact, the problem which triggered the work presented 
in this paper was the problem of reconstruction of straight line image tracks in pixelized single-photon detectors 
with rectangular pixels, designed for astro-particle physics experiments . Similar problems are often encountered in 
high-energy particle physics, when dealing with the reconstruction of particle tracks, and in the field of imaging by 
means of pixelized detectors. In astrophysics color-color diagrams typically lead to similar problems (see section V A|. 
The method proposed in this paper has been also applied irP. Moreover the SWM is often assumed whenever the 
errors of the data points are unknown and there is no reason to assume that the errors in one variable are negligible 
with respect to the errors in the other one. In this case the common unknown value of the error can be estimated by 
the fit. 

The LSFSL-SWM problem can be always reduced, by a suitable rescaling of the variables, to an equivalent problem 
with the new dimensionless variables having equal error^. The sum of the distances between the measured data 
points and the best-fit straight line can thus be minimized leading to a purely geometrical problem. It is possible to 
solve this problem in an exact way and the results are very simple. 

The purpose if this paper is to present a complete set of explicit results for the LSFSL-SWM. The results from 
the physics literature known to the author are quickly reviewed, and other new results are derived, including very 
simple analytic formulas for the standard errors of the parameters. It is shown that parameterizing the straight line 
with the angle with respect to one of the axes plus a second dimensional parameter leads to very simple results, 
having simple transformation properties under roto-translation of the Cartesian Coordinate System. It is shown that 
variances and covariance of the parameters can always be expressed in terms of the standard error of the angle plus 
purely geometrical quantities. All the results and the accuracy of the expressions for the standard errors have been 
cross-checked with extensive Monte Carlo simulations. The results are also compared with the results of the ordinary 
least squares (OLS) fit, the case of significant errors in only one of the variables, which will be called, in shorlP, 
OLS-y:x fit (error on the y variable only) and OLS-x:y fit (error on the x variable only). A simple criterion for 
neglecting the error in one of the two variables is derived. 

Emphasis of this paper is on the analytical, computational and practical aspects of the problem, not on its rigorous 
statistical treatment. All the formulas needed for a real implementation in a programming code are derived and 
discussed. 

The outline of the paper is as follows. In section |ll] the problem and hypotheses are presented. Results readily 
available in the physics literature are summarized in section |III| All the new results are presented in Section |IV| 
Section|V]presents a few case studies. Section [VT| summarizes some of the results obtained by Monte Carlo simulations 
carried on both to cross-check the formulas and to evaluate the accuracy of the standard errors. Finally the 

appendixes collect all the calculations. 



II. HYPOTHESES AND THE SWM 



For the sake of simplicity let us introduce the symbol for the average of the generic random quantity q, (q), the 
symbol a[q] for its standard deviation and the symbol ppq for the correlation coefficient of the two random variables 
p and q. 

Consider a set of N measured data points, described by the x and y coordinates of a suitable Cartesian Coordinate 
System: {xk;yk} , {k = 1,...,N). Assume that the measured data points follow the SWM, that is: have equal 
standard errors, separately in both variables, ax = clxk] and ay = a[yk], and that the errors on the x and y variables, 
for each measured data point, are uncorrelated: 

{xk ± as; yk±ay} (fc 1, . . . , N) p^y = . (1) 

It is assumed that each measured data point is a random sampling from a random distribution associated to each 
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true data point, {X^; Yfe}, with true values lying on a unknown straight line, to be determined: 

{Xk^Yk} AXk + BYk+C = Q (fc = l,...,7V) (2) 
{xfc = Xfc + ejr; yfe = Yfc + ey} with (e^) = cr[e2] = ctj; (cy) = a[ey\ = (Jy , (3) 

where and ejj are random variables, often, but not always, Gaussian random variables. 

The problem can be always reduced to an equivalent problem^ with identical errors in both variables, by rescal- 
ing, via the standard errors, to dimensionless variables, and multiplying, for the sake of generality, by a common 
dimensionless factor, r: 

Xk — > Xk = T (xfe/cr^) => (7[xk] =T yk — > yk = T ivk/cTy) =^ a[yk] = t . (4) 

The original variables, Xk and yk, will be called the raw variables, as opposed to the Xk and yk, the (re-scaled) variables. 
The above transformation will be always silently assumed in the rest of this paper. In every physics problem this 
transformation is just a change of the units of measure of the two variables, using and ay as the new units of 
measure, leading to dimensionless variables. One might choose t — 1, but it is preferred to leave a generic r in order 
to have a better understanding of the final formulas and deal with the case of a common but unknown error. 

After the above transformation, the sum of distances between all the measured data points and the straight line, 
(the error function) can be minimized, as a purely geometrical problem. Sector the precise probabilistic/statistical 
discussion. 

Note that there is no point in discussing any ambiguity of the best-fit line under change of scale of the coordinates. 
In fact only when the errors in both variables are equal one is allowed to minimize the distance between the mea sured 
data points and the straight line, otherwise one needs to take into account the error ellipse, see for instanc^^. 

The OLS-y:x/OLS-x:y fit can be recovered by letting crj — > 0/ay — >■ in the final results. 

Finally, let us introduce the following notations for the variance and covariancc of the set of N measured data points 
as whole, defined without the Bessel correction factor N/ {N — 1), as it is useful in the formulation of the LSFSL 
problem: 

V, = {x^)~{xf Vy = {y')-{yf a,y = (xy) - (x) (y) AV = V, - Vy . (5) 

The above quantities, Vx, Vy and C^y, refer to the set of measured data points as a whole, describing their spatial 
distribution in the xy plane: they are obviously distinct from the spread of one single measurement around its true 
value. 

In order to avoid any confusion, the symbols Vx, Vy and Cxy, are used in this paper to refer to the set of measured 
data points as a whole, while the symbols cr^, ct^ and Pxy = are used to refer to the variances and correlation 
coefficient of the two-dimensional random variable associated to each specific measured data point. 

The relations between variances and covariance in the raw and re-scaled coordinates, necessary for taking the 
no-error limits, are obviously: 

Vx^T^^ Vy^r^^ Cxy^r'^ . (6) 



III. THE SLOPE/INTERCEPT PARAMETRIZATION 

The least squares fit to the straight lines 

y PyX + Qy py = tan[6':r] or x = pxy + qx Px = ta.n[9y] , (7) 

in terms of the slopes, Py/px, intercepts, qy/qx, and angles 9x/0y, with respect to the x/y axes, is appropriate in many 
physics problems, for instance when fitting a position as a function of time, so that the slope (velocity) can be zero 
but not infinity. 

Some of the most well-known textbook^SHUl ygg this parametrization to describe straight lines and present the 
so-called effective variance methoJ^, to deal with measured data points having significant errors in both variables. 
One exception is textboolP, which presents the exact solution within the SWM. 

The dimensionless error function to minimize, the sum of distances between all the measured data points and the 
straight line, is 

2r 1 1 (PyXk+qy-Vkf 2r 1 ^ {PxVk + qx " Xkf 

Xv[Py,<ly]= -^2^ 7—^ or Xx[Px,qx] = — ^ , (8) 

k=i k=i '^^i'x 
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which will be called x^; regardless of its statistical properties. 

Minimization of the error functions leads, after some mathematics, to the expressions for the slopes and inter- 
cepts which can be found in some of the literatur^S^El, xhe solution, in both the y versus x and the x versus y 
representations, is: 

for C^y ^ with a = Sign[C,^] = ^ (9) 



Ay = ' " Py = Ay + aJl + Al^— ' : qy ^ {y) - pyi^x) (10) 



^^xy ^^xy 

Vx - Vy , (K - Vy) + j^Cly + {l^Vf 

Ax^^-^ Px^Ax + a^lVM= = {x) - Px{y) (H) 

Z^xy Z^xy 

with PyPx = 1 and 9x + dy — 7r/2 . (12) 
Note that the sign in front of the square-root is unambiguously determined by a. 

The above relation PyPx = 1 is between the two different representations of the same best-fit straight line. It should 
not be confused with the relation between the Py and Px slopes as determined by the OLS-y:x/OLS-x:y fit, giving the 
correlation coefficient, r, of the set of data-points as whole (see for instanc^^^: PyPx = r^- 

Whenever Cxy = the fitted line is parallel to either the x oi y axis. The sign of the slopes is always the same 
as the sign of the Cxy The straight line perpendicular to the minimum-distance straight line maximizes the error 
function; its slope can be found by putting, in the above formulas, a minus sign, instead of a plus sign, in front of the 
square root. 

A simple trigonometric transformation allows to transform equations [TO] and [TT] into an even simpler equation for 
the angles 9x and dy: 

ta.n[29x] - ^ for 7^ or tan[26';,] = - ^ for A^, 7^ . (13) 

Ay Ax 

It is useful to re-write the above equations [TOl and [TT] as a function of the raw variables, which also allows to take 
the hmit of the OLS-y:x/OLS-x:y fit. One find&^> 

p,.tan[.,]= J(A,+.,;rT:4^ , (14) 

'^IVx - aWy 



Px ^ tan[%] = ^ + ayr+Af) Ax ^ ! ^ „^ ^ . (15) 
Oy \ J zaxiyyUxy 

The standard errors on the slope and intercept will be derived in section |IV G| 



IV. THE ANGLE/SIGNED-DISTANCE PARAMETRIZATION 
A. Parametrization of a straight line in the plane 

In many problems the x and y coordinates are equivalent, such as, for instance, in straight line fitting of pixelized 
images with rectangular/square pixels. Especially in these case a parametrization of the straight line different from the 
most common slope/intercept parametrization may be useful, as the direction of the straight line is better identified 
by the angle with respect to one specific direction in the plane, for instance the x axis, than by the slope. In fact 
both the slope and intercept may tend to infinity for straight lines passing near the origin and nearly parallel to the 
y OT X axis (depending on the representation used from equation [7| . 

The details of the geometrical aspects of the parametrization of a straight line in the plane which are useful for the 
understanding of the results are described in this section. 

Any straight line in the xy plane passing by the point R = {X; Y} can be represented in terms of the direction 
unit vector, n = {cos[9]; sm[9]}, by the parametric equation 

r[t]^R + nt n^{cos[0];sin[0]} R^{X;Y} ft! ^ ^ ' (!«) 

\y[t\ = y + tsm[6'J 

Two real parameters identifying in a unique way all straight lines in the plane can be chosen as: 
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• the angle, 9, between the straight hne and the x axis: 

— 7r/2 < < +tt/2 cos[6'] > the range of 9 is just one possible conventional choice ; (17) 

• the signed-distance c, given by the z component of the vector product between n and r, whose absolute value 
gives the distance between the straight line and the origin: 

— oo < c < +0O c = Bz ■ {n X R) for any point R belonging to the straight line ; (18) 

in fact the magnitude and sign of c is also given by the vector product between n and the vector position of any 
point on the line, R; the sign of c is always the same as the sign of Qy. 

The geometry is represented in Figure [T] 




FIG. 1: The angle/signed-distance parametrization. 

Converting into the Cartesian representation one finds: 

xsin[6'] - ?/cos[6'] + c = (19) 
c = —X sin[9] + Y cos[9] for any point R = {X; Y} belonging to the straight line (20) 
{x - X) sm[9] - (y-Y) cos[9] = . (21) 

Of course, other parameterizations are possible, but the above parametrization has the virtue of describing bi- 
univocally all straight lines in the plane in terms of two real non-singular parameters with a clear and direct geometrical 
interpretation. Moreover the two parameters 9 and c have simple transformation properties under translation and 
rotations in the xy plane, at variance with the slope and intercept. 



B. Determination of the best-fit line 

The error function to minimize using the angle/signed-distance parametrization is: 

1 ^ 

X^[9, c]= (^'^ ^i'^t^] - + '^)' ' (22) 

k=l 

which will be called regardless of its statistical properties. 

Note that the problem is a non-linear least squares problem, so that the nice general properties of linear least 
squares problems (see for instance^ are not guaranteed. In particular the exact contour of the confidence region in 
the parameter space is not ellipticaP^. 

All the details of the minimization of the error function in equation [22] are summarized in appendix |X] where all 
the following results are derived. 
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Equating to zero the derivative with respect to c of equation 22 in order to look for stationary points of the error 
function, immediately provides the condition that any straight line leading to a stationary point of the error function 
passes by the centroid of the data points, C = {{x); (y)}, (as it is the case for the OLS-y:x/OLS-x:y fit): 

i^^^^M^O— >c = -(2;)sin[^] + (t/)cos[^] . (23) 
oc 

The error function in equation |22| becomes a function of 6 only, after using equation[23]to replace the signed-distance. 
It can be written as: 



N = ^E((^'=~<^))^i'^[^]-(yfc-<y))^°"[^])' = "^-^i'^'[^H"^/Cos2[0]-2C^^ . (24) 

fc=i 

The values of 6 giving the stationary points are (see appendix |A] for the derivation of the results) : : 

AF = =^ cos[26i] = =^ ^"I+^i "^^^ 
AV^O =^ tan[20]=^ 0~ ^ 1 arctan[ ^] + g ^ q€Z (26) 



minimum: choose 9 such that 16*1 < tt/2 with the same sign as C; 



xy 



In the often encountered case that the common error on the data points, t, is unknown, it can be estimated from 
the data. It can be showrflS'that, in the linear least squares case, the residual sum of squares at the minimum, divided 
by the number of degrees of freedom, N — 2 (as two parameters are determined by the minimization), is an unbiased 
estimate of the unknown error r: 



q2 ^min /'07^ 



C. Formulas for the standard errors 



The simple exact analytic formulas for the standard errors and covariance of the fit parameters are presented and 
discussed in this section. Ex act means tha t they can be derived without any approximation from the standard formula 
for the propagation of error^SHEnEHHIIZl^ which is derived using the first-order truncated Taylor series expansion of 
the function and taking expectations values. Approximate formulas were published iiP. 

A long and tedious, but straightforward, direct calculation, applying the standard formula for the propagation of 
errors to equation 13 leads to the following results for the standard errors (see the appendix |A] for details of the 
calculation) : : 



n = {cos[6l];sin[6']} C={{x);{y)} Z = (x) cos[d] + {y) sm[§] = n ■ C 



Varfei 



N [AVY + ACly 

2 



Var[c] = {5cY = — + [Z 56)' 
CoY[e,c] = -Z{6ef . 



(28) 
(29) 

(30) 
(31) 



Equations 
in terms of t 



29 30 and 



le standarc 



31 generalize the well-known results for the OLS-y:x/OLS-x:y fit. They can all be expressed 



error on the angle plus the geometrical term Z . 
Z is interpreted as the signed-projection of the position vector of the centroid, C, onto the direction unit vector, n, 
of the best-fit line. The absolute value of Z is the distance between the centroid and the straight line perpendicular 
to the best-fit line and passing by the origin. 

All the above standard errors, equations |29[ |30| and |31[ must be invariant under rotations in the xy plane as c is 
invariant while 9 is just shifted under rotations. However this property is not obvious from equations [29] [30] and [31] 
it is demonstrated in section ITV El 

The standard error in equation [29] must be also invariant under translations in the xy plane, because 6 is invariant 
under translations, while c is not. This property is obvious from equations |29[ [30] and [3Tj as only Z is not invariant 



under translations in equations 29 30 and 31 
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Equation 30 shows that the error on the signed-distance c is minimized if the origin of the Cartesian Coordinate 
System is set at the centroid of the measured data points, so that Z = (the same property apphes in OLS-y:x/OLS-x:y 
fit). In this case the error on the signed-distance c is just the expected t/^/N, that is the error on the determination 



of the mean values (x) and (y), according to equation 23 Moreover if Z = the estimates of and c are uncorrelated. 
In the general case equation [30] is the sum in quadrature of a term coming from the uncertainty in the position of the 
centroid, the first one, plus a term coming from the uncertainty in the angle, amplified by the distance \Z\. 

It should be emphasized, as pointed out in some of the literature, that the formulas for the standard errors should 
be evaluated, in principle, via the true data points and not via the measured data points, as the standard error 
expressions are derived from a Taylor series expansion about the true data points. In practice the measured data 
points are normally used to evaluate the standard error expressions. The accuracy of this approximation is studied 
iiP, via Monte Carlo simulations. In principle, after the fit, one might want to correct the measured data points to 
estimate the true data points values, improving the calculation of the standard errors. However it is shown in sectiorPl 
that this is normally not necessary. 

D. Bias 

As the problem is a non-linear least squares problem, un-biasedness of the estimates is not guaranteed. 

Bias of the estimated parameters, if any, is invariant under roto-translations of the Cartesian Coordinate System. 

It is easy to show that the estimator for the angle 9 is un-biased. Imagine to set the origin of the Cartesian 
Coordinate System at the unknown position of the centroid of the true data points, with one axis along the true 
straight line. It is then clear, by symmetry, that the resulting probability distributions of all the measured data points 
are symmetrical around the unknown true straight line. Therefore also the distribution of the fitted angle will be 
symmetrical around the true value and the estimation of 9 is un-biased. 

Similarly, the estimator for the signed-distance c is un-biased as well. 

E. The covariance matrix and Principal Components Analysis 

Principal Components Analysis (PCA) provides a better insight into the results of the LSFSL-SWM and help to 
derive some results in a simple way. It relies on orthogonal transformations of random variable&^°. PCA and its 
relation with LSFSL-SWM, discovered by K. Pearsorfi^, are extensively discussed ini2l. 

In the rotated Cartesian Coordinate System (x; y) — (^; rf), such that 9 is the angle made by the ^ axis with the x 
axis, one has: 

^ = +j: cos[e'] + y sm[9] rj ^ -x su\[9] + y cos[6'] , (32) 
and the variances of the set of data points are readily calculated: 

[V^^V^ cos2 [9] + Vy sin2 [9] + 2C^y sin[0] cos[0] 

I F,, = K sin^ [9] + Vy cos2 [9] - 2C^y sin[9] cos [9]=V^ + Vy-V^ . ^ ' 

Therefore requiring that the variance, V^, of the set of measured data points as a whole is minimized, leads to 
exactly the same results as minimizing equation |24[ In fact the LSFSL-SWM minimizes the sum of the distances 
between the measured data points and the straight line at angle 9 from the x ax:is, namely the ^ axis. 



Equations 33 show the well-known result that the sum of the variances is invariant under rotations, as it is the trace 
of the covariance matrix. Therefore minimizing Vr, is equivalent to maximizing V^. Som e standard linear algebra, see 
for instanccJi^, shows that the direction ^ maximizing the variance in equation 33 can be found by diagonalizing 



the covariance matrix: 



(34) 



Therefore diagonalizing the covariance matrix leads to maximize the variance along the ^ axis and minimize the 
variance along the rj axis: the largest eigenvalue, A-|„, is the variance along the ^ axis and the smallest eigenvalue, 
A_, is the variance along the 77 axis. The two orthonormal eigenvectors give: the direction maximizing the variance, 
namely the ^ axis, (the eigenvector corresponding to the largest eigenvalue, {cos[0]; sin[0]}); the direction minimizing 
the variance, namely the rj axis, (the eigenvector corresponding to the smallest eigenvalue), corresponding to minimize 
the error function in equation |22| 
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The following results can be obtained: 



X±^ -[{V^ + Vy)±^ (Ay)' + ACl^ ) (35) 



0<A_<A+ \+ + \^=V^ + Vy \+ - \- = {/SVy + 'iCly . (36) 
The case of perfect linear correlation corresponds to C^y — VxVy, that is zero determinant of the covariance matrix. 

F. Some properties of the standard errors 



The discussion in section IV E gives some insight into the understanding of the properties of the standard errors in 
equations [29l |30] and |3l] 

The expressions for the eigenvalues of the covariance matrix, equations |35[ show that the standard errors in equa- 
tions [291 [30] and [31] are invariant under rotation in the xy plane, as they can be written in terms of the eigenvalues of 



the covariance matrix and the signed-distance Z. 

In the rotated Cartesian Coordinate System (C; 77) , one has: 

= A+ K = A_ [Sef = — ^ ^ . (37) 

In the often encountered case that the measured data points are highly linearly correlated, the covariance can be 



approximately written in terms of the variances, C^y ~ V^Vy, and equation 29 simplifies to: 

2 t2 1 t2 1 

= Iv viTk; » ■ 

The above result for highly linearly correlated measured data points has a very simple interpretation. It shows that 
the standard error on the angle can be expressed as the ratio between the single data point standard error, r, and 
the square-root of the variance along the best-fit straight line of the set of measured data points as a whole, \/V^, 
divided by the square-root of the number of measured data points. Therefore the standard error on the angle can be 
interpreted as the transverse size of the single measured data point, r, divided by an effective length of the measured 
data points set, given by ^/V^, similarly to the definition of the radian, all divided by the square-root of the number 
of measured data points. 

This result can be used to easily estimate the error expected from the LSFSL-SWM, for highly linearly correlated 
data. Moreover it makes quantitative the naive expectation that the data points at the two extremes of the straight 
line have more importance for the fit, as they increase the variance decreasing the error on the angle. In fact: 
SV^ = 2{C,~{0)6^,/N. 

G. Standard errors on the slope and intercept 

Using the results derived in this section it is possible to find the formulas for the standard errors on the slope and 
intercept: 

{5pf^{l+p')\6ef (39) 



{sqf = (1 ( ^ + - ^vf) = (1 ( ^ + {-r (i +?')) 



(40) 



See appendix [C] for the details. 



H. Comparison with the OLS-y:x/OLS-x:y fit 



The invariance under rotation of the error function (as it is a sum of distances) can be used to compare the results 
with the OLS-y:x/OLS-x:y fit. 
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The error functions in equation [8] clearly show that the results of the fit will tend to the results for the 
OLS-y:x/OLS-x:y fit whenever py ~ 0/px — 0, as in this case, the distance is measured in a direction which is 
both perpendicular to the straight line and parallel to the y/x axis. 

In fact, after determining the best-fit straight line one can apply a roto-translation to bring the origin of a new 
Cartesian Coordinate System, {^;?7}, coincident with the centroid of the measured data points, and the ^ axis along 
the best-fit line. In this new Cartesian Coordinate System the error function in equation 22 is exactly the same as it 
would be for the OLS-ry:^ fit. The covariance of the set of measured data points, C^,,, is zero and the best-fit slope 
and intercept are both zero, in both cases by construction of the Cartesian Coordinate System l^;??}. 

A more detailed discussion can be found in appendix [B] 



I. OLS - negligible error in one of the variables 

In the case of the highly linearly correlated measured data points of equation [38j it is easy to find a simple criterion 
to assess whenever the errors in one of the two variables can be neglected, the OLS fit. 
After re- writing equation [38] in terms of the raw variables one finds: 



1 



1 



1 



N V^ + Vy 



N Vy 



(42) 



Note, however, that equation 42 does not give the limiting standard error on the raw angle 9, which is the interesting 
quantity, but the one on the re-scaled angle 9. 

shows that the error on x/y is negligible depending on the relative 



The complete analysis, in the appendix B 



values of Vx/cr- versus Vy/a~ 



X y 



OLS-y:x Spy 
OLS-x:y Sp^ 



1 (Ty 



(43) 



N ^Vy 



V. CASE STUDIES 



A. B. C. Reed's example 

Consider B. C. Reed's example pSl^ which was analyzed irP by using approximate expressions for the standard 
errors. These data derive from a real, physical situation: calibrating the colors of globular star clusters as a function 
of their spectral types. The abscissa represents the difference between the ultraviolet and yellow light magnitudes of 
the clusters and the ordinate represents the difference between the yellow and infrared magnitudes. 

The results are summarized in Figure [2] and Table [l] . The results for the LSFSL-SWM are in perfect agreement 
with the paper of the author while the approximate treatment of papeil^l leads to slightly reduced error estimates'^ . 



OLS-y;x 


Pv = 


f0.9311± 0.1299 


Qy = 


-0.1501 ±0.1192 


OLS-x:y 


Py = 


fl.3839± 0.1930 


ly = 


-0.5631 ±0.1768 


LSFSL-SWM 


Py = 


fl.l668± 0.1704 


qy = 


-0.3652 ±0.1561 



TABLE I: Best fit results (slope, intercept): OLS-y:x, LSFSL-SWM and OLS-x:y. 



Since the data are given with only two significant figures, the number of digits quoted is not justified. However, 
following B. C. ReecP^, the point is to provide figures against which others can compare the results of their own 
algorithms as if the original data are regarded as exact. 

The results of PCA are shown in Figure [3] 
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FIG. 2: Data points and the three best fit lines: OLS-y:x (the smallest slope), LSFSL-SWM (the intermediate slope) 

and OLS-x:y (the largest slope). 




FIG. 3: Data points and results of PGA: the ellipse has its semi-axes along the directions of the eigenvectors and the 

semi- axes have the same length as the eigenvalues. 



B. A light-detector 



Consider a CCD-like light detector, where the image produced by a suitable optics is focused (a camera, for 
instance). Suppose it must be used to identify straight light tracks, produced by a moving point-like light source in 
the night sky, for instance to observe meteors, airplanes or extensive air-showers produced by cosmic radiatiorP. Let 
the typical length of the track be ctl , defined as the standard deviation of the data points along the track. 

Assume, firstly, that the combination of the intrinsic spread of the point-like light source and the point-spread 
function of the optics gives, on the photo-detector plane, a standard deviation — Oy = t much larger than the pixel 
size, so that binning effects can be neglected. The expected angular resolution can be estimated as (equation 38 ) : 



T 



T 1 



(44) 



As a second example consider the case that the pixel size, d, is much larger than the transverse size of the image track 
on the photo-detector plane. The uncertainty on the measurements can be taken as r = d/-\/l2, assuming a uniform 
probability distribution inside the pixel. Assuming uncorrelated measurements, the expected angular resolution can 
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be estimated as (equation 38 ) : 



66 



T 

N 



1 



1 



N V12 <JL 



(45) 



Both analytical relations can be used to optimize the photo-detector parameters, given the desired angular resolu- 
tion. 



C. Kinematics of a image traclc in a photo-detector 

The explicit analytical formula 29 was used, in referenc^, in the design of a photo-detector to observe the image 
track of the extensive air-showers produced in the atmosphere by ultra high energy cosmic radiation. The explicit 
analytical formulas, written in terms of the photo-detector parameters, were then used to optimize the photo-detector 
design. 



VI. TOY MONTE CARLO SIMULATIONS 
A. Setup of the toy Monte Carlo simulations 



In order to cross-check and evaluate the accuracy of the formulas for the standard errors derived in this paper, 
extensive Monte Carlo simulations' were carried on as follows. 

1. Straight lines were simulated in the xy plane, with uniformly distributed random angle 6 and a Gaussian 
distribution of the signed-distance, c, with zero mean and a standard deviation oi D = 1. For each straight line 
a fixed number, N, of true data points was simulated, on the straight line, with a uniform random distribution 
along a segment of length L = 1 of the straight line. 

2. Each true data point was converted into a simulated measured data point (i.e. simulated measurements) 
according to a bi-dimensional Gaussian distribution centered at the true data point and having equal fixed 
standard errors, t = (7^ — <^y, and zero correlation coefficient. The best-fit straight line was then calculated. 

3. For each simulated straight line and its set of true data points, the procedure at step 2) was repeated for a 
number of times; each repetition is called one iteration and, typically, the number of iterations is nj sa 1000. 
Statistics was then collected of the fitted angle/signed-distance results, for one specific true straight line and its 
set of true data points. 

4. The procedure at steps 1), 2) and 3) was repeated by simulating a number of different straight lines; each 
simulation of one straight line is called one run and, typically, the number of runs is nj^ w 1000. For all the 
runs, statistics was collected of the fitted parameters with respect to the known true parameters of the simulated 
straight line and its set of true data points. 

5. The procedure at steps 1), 2), 3) and 4) was then repeated for different values of the number of measured 
data points, N, and of the standard error of the single data point, r. The different values for N and r where 
chosen in such a way to keep an approximately constant ratio between one value and the next one, inside a 
pre-defined fixed range. Therefore different sets of parameters, {N,t), were simulated, with: 3 < < 100 and 
0.001 < r < 1, as follows: 



N = {3; 4; 5; 6; 8; 10; 15; 20; 30; 50; 100} ; (46) 

^ j coarse-scan in r: r = {1; 0.1; 0.01; 0.001} cx lO"'' 

fine-scan in r: r = {0.1; 0.0631; 0.0398; 0.0251; 0.0158; 0.01} oclO-*^/^ 

The problem is invariant under a common rescaling of the two coordinates, so that the results are expected to 
depend only on the ratio t/L. 

The largest simulated value of the data point error, t/L = 1, would give so large errors with respect to the length 
of the segment of the straight line that it is a non realistic case. It has been simulated to verify that, in this case, the 
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approximation of the standard formula for the propagation of errors fails and to determine under what conditions the 
formula becomes an accurate estimate of the standard error. 

It would be also possible to simulate straight lines with absolute values of the signed distance, c> L = 1. However 
the standard error formulas, equations [29] [30] and [31] show that it is better to set the origin of the Cartesian Coordinate 
System as close as possible to the centroid of the measured data points, in order to have c ~ 0, so that Z ~ and 
both the error on the signed distance and the covariance between the angle and the signed-distance are minimized. 
This is of course always possible, by means of a suitable translation, and therefore there is no need to simulate larger 
absolute values of the signed-distance. 

For both the angle 6 and the signed-distance c, the three following quantities were calculated. 

1. The standard deviation of the distributions of the fitted angle and signed-distance, a[0] and (t[c], were calculated 
for every run; for any set of parameters the average was taken over all the runs: {^[0]) and ((7[c]). 
This is what is usually defined as the standard deviation of the estimator. 



Equations 29 and 30 were used, for all iterations of every run, to calculate the standard errors from the simulated 
measured data points (i.e. simulated measurements) and the median was taken over all iterations of every run, 
M[d9] and M[(5c]; for any set of parameters the average was taken over all the runs: {A'I[S9]) and {M[Sc]). 

This is an estimate of the typical error one would calculate from a real set of data points. The median is used 
in order to provide a robust estimation, as outlying large values may show-up for certain combinations of true 
data-points (for instance all true data-points clustered one close to the other). 



3. For each run the standard errors were computed from the true data points via equations [29] for 69o, and 30 for 
(5co; for any set of parameters the average was taken over all the runs: {56o) and {5co). 

For any given straight line and any set of true data-points on it, 69o and Scq are the expected standard errors 
calculated by the standard formula for the propagation of errors. Therefore {69o) and (Scq) are a sort of reference 
values, setting the expected magnitude of the standard errors, for any given set of parameters, and useful as 
normalization factors to compare different set of parameters. 

Some of the results obtained by the extensive Monte Carlo simulations are summarized in the rest of this section. 



B. Results of the toy Monte Carlo simulations 



The results of the different simulations were re- normalized by r, in order to compare results for different values of 
T, and are shown as a function of N for different values of t/L. 



1. Results for {59o) and (Scq) 

All the simulations for {69q) and (Scq), after normalizing for r, show a very similar behavior as a function of N, for 
the different values of t/L, as shown in Figure [4] 

The behavior as a function of N is fitted by ^ l/\/N + 4, for {S9o), and by ^ ^/\/N -f 1, for (Jcq) (the result of 
the fit not shown in the Figures). 
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FIG. 4: Results of the simulations for {S9q) and (Scq), as a function of TV, for t/L = {1; 0.1; 0.01; 0.001}. The four 
curves on each plot are well superimposed and barely distinguishable. Errors bars are very small and barely visible 

as well. 
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2. Results for {a[0]) and {a[c\) 



All the simulations for ((5'[0]) and (^[c]), after normalizing for r, show a very similar behavior as a function of N, 
for the different values of t/L < 0.05, as shown individually in Figure [5] and all together in Figure |6j 

Moreover, for values t/L < 0.05, the behavior as a function of N is fitted by ~ + 4, for {(j[0]), and by 

~ + 1, for ((t[c]) (the result of the fit not shown in the Figures). 

In order to study the behavior for 0.01 < t/L < 0.1, a finer scan has been done in this interval; the plots 
are individually shown in Figure [?) for t/L = {0.0631; 0.0398; 0.0251; 0.0158}, and all together in Figure [sj for 
t/L = {0.1; 0.0631; 0.0398; 0.0251; 0.0158; 0.01}. 

Afterward, in order to get rid of the variability associated with the random straight line and random set of true 
data-points on the straight line, the values of {(j[9]) and {a[c]) were normalized to (^^o) and (Scq). 

For r < 0.02, no statistically significant difference with {S6q) and (Sco) has been found and the ratios {a[9]) / {S9q) 
and ((t[c])/((5co) are compatible with one, as shown in Figure |9] 

On the other hand the behavior of the ratios {a[9]) / (SOq) and {d-[c]) / {6co) for values oi t/L > 0.02, as shown in 
Figure 10 shows a significant departure from the expected {d9o) and (Scq), but not larger than « 5% for r < 0.1. 

As a general trend, the ratios tend to one as t/L decreases, as expected thanks to the improvement of the ap- 
proximations made to derive the standard formula for the propagation of errors. Moreover the ratios tend to one at 
small TV anyway because at small A'^ the values of the numerator and denominator became large with respect to their 
difference. 
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FIG. 5: Results of the simulations for {a[9]) and ((t[c]), as a function of N, for t/L = {1; 0.1; 0.01; 0.001}. Errors 

bars are very small and barely visible. 
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FIG. 6: Results of the simulations for {a[9]) and {&[£]), as a function of A^, for t/L = {1; 0.1; 0.01; 0.001}. The two 
almost horizontal curves on each plot correspond to t/ L — 1. The two curves with the smallest t/L are well 
superimposed and barely distinguishable. Errors bars are very small and barely visible. 



17 




FIG. 7: Results of the simulations for {a\9]) and (^[c]), as a function of N , for 
t/L = {0.0631; 0.0398; 0.0251; 0.0158}. Errors bars are very small and barely visible. 
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FIG. 9: Results of the simulations for the ratios {a[d]) / {d9o) and {d-[c]) / {6cq), as a function of N, for t/L = 0.0158. 
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FIG. 10: Results of the simulations for the ratios {a[9]) / {S9o) and {a[c]) / (Scq) , as a function of N, for t/L — 0.1. 
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3. Results for {M[5e]) and {M[5c]) 



All the simulations for {M[66]) and {AI[Sc]), after normalizing for r, show a very similar behavior as a function of 
iV, for the different values of t/L < 0.05, as shown individually in Figure 



11 



and all together in Figure 12 



4, for {M[50]), and by 
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for 



Moreover, for values t/L < 0.05, the behavior as a function of N is fitted by 
~ + 1, for (M[(5c]) (the result of the fit not shown in the Figures). 

In order to study the behavior for 0.01 < t/L < 0.1, a finer scan has been done in this interval; the plots 
are individually shown in Figure 13 for t/L = {0.0631; 0.0398; 0.0251; 0.0158}, and all together in Figure 
t/L = {0.1; 0.0631; 0.0398; 0.0251; 0.0158; 0.01}. 

Afterward, in order to get rid of the variability associated with the random straight line and random set of true 
data-points on the straight line, the values of {M[59]) and {M[5c\) were normalized to {59q) and ((5co). 

For T < 0.02, no statistically significant difference with {60q) and {Scq) has been found and the ratios {M[5ff\) / {60q) 
and {M[5c]) / {5ca) 
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are compatible with one, as shown in Figure 
On the other hand the behavior of the ratios {M\59]) / {59{)) and {M[5c\) / {5cq) for values of t/L > 0.02, as shown 
in Figure 16 shows a significant departure from the expected {S9q) and {Scq), but not larger than « 10% for r < 0.1. 

As a general trend, the ratios tend to one as r/L decreases, as expected thanks to the improvement of the ap- 
proximations made to derive the standard formula for the propagation of errors. Moreover the ratios tend to one at 
small N anyway because at small N the values of the numerator and denominator became large with respect to their 
difference. 
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FIG. 11: Results of the simulations for {M[6e]) and {M[Sc]), as a function of TV, for t/L = {1; 0.1; 0.01; 0.001}. 

Errors bars are very small and barely visible. 
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FIG. 12: Results of the simulations for {M[5d]) and {M[Sc]), as a function of iV, for r/L = {1; 0.1; 0.01; 0.001}. The 
two almost horizontal curves on each plot correspond to r/L = 1. The two curves with the smallest t/L are well 
superimposed and barely distinguishable. Errors bars are very small and barely visible. 
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FIG. 13: Results of the simulations for {M[59]) and {M[Sc]), as a function of N, for 
t/L = {0.0631; 0.0398; 0.0251; 0.0158}. 
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FIG. 14: Results of the simulations for {M[S9]) and {M[Sc]), as a function of N, for 
t/L = {0.1; 0.0631; 0.0398; 0.0251; 0.0158; 0.01}. Errors bars are very small and barely visible. 
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FIG. 15: Results of the simulations for the ratios {M[S9]) / {S9q) and {M[5c]) / (Scq), as a function of N, for 

t/L ^ 0.0158. 
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C. Conclusions from the Monte Carlo simulations 



The main results of the Monte Carlo simulations presented in the previous sections can be summarized as follows 

and 



• The standard errors, equations 29 and 30 calculated from the true data points, ((5(.)q), scale as ^ l/\/N + A, 
for 6, and as ^ + 1, for c, at fixed t/L. Normalization by r provides, for t/L < 0.02, a universal curve 

as a function of N. 



• Both the standard deviation of the estimators of the parameters, 9 and c, ((t[(.)]), and the standard errors, 
calculated from the measured data points, {M[6{.)]), as long as t/L < 0.05, scale, at fixed 



equations 
t/L, as ^ 



29 



1/ 



and 

N 



30 



4, for 9, and as ^ l/-\/iV + 1, for c. Normalization by r provides, for t/L < 0.02, a universal 



curve as a function of N. 



The standard errors, equations 29 and 30 calculated from the measured data points, {M[S9]) and {M[6c]), are 
reliable estimates of the standard deviation of the estimators of 9 and c, {&[{.)]), to within sa 10% for t/L < 0.1. 
In particular, to within « 10%, it shall not be necessary to re-evaluate the standard errors using the estimated 
true data points values after the best-fit straight line has been determined. 



For t/L < 0.02, the standard errors, equations 29 and 30 both calculated from the true data points, {5{ 



>0h 

and calculated from the measured data points, (M[(5(.)]), and the standard deviation of the estimators of the 
parameters, ((t[(.)]), do not show any difference, within the statistical uncertainty of the Monte Carlo simulations. 

The following additional conclusions were obtained by the Monte Carlo simulations, and stated without showing 
explicit evidence in this paper. 



• Lack of bias of the two estimators, to within the statistical uncertainty. 

• Excellent normality of the distribution of the 9 estimator, for all the simulated parameters, to within the statis- 
tical uncertainty, according to common statistical tests, such as Cramer-von Mises and Kolmogorov-Smirno\^. 

• Excellent normality of the distribution of the c estimator, whenever t/L < 0.01, to within the statistical 
uncertainty. For larger values oi t/L strong deviations from normality start to be very significant for small 
values of N and less significant for large values of iV; for instance, Cramer-von Mises and Kolmogorov-Smirnov 
tests give a P-value less than « 0.001 whenever t/L — 0.1 and N < 37. The distributions, deviating from a 
Gaussian shape, show a positive excess kurtosis. 



VII. CONCLUSIONS 



Simple formulas for both the best-fit parameters and their standard errors in the LSFSL-SWM have been derived, 
using the angle/signed-distance parametrization of the straight line, and validated with Monte Carlo simulations. 

Several properties of the standard errors were derived and investigated. The standard errors in the slope/intercept 
parameterization were derived. A simple relation for the case of highly-correlated measurements was derived. The ex- 
tent to which errors in one of the variables can be neglected is quantified and the limiting case of the OLS-y:x/OLS-x:y 
fit was studied. 
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Appendix A: Derivation of the best-fit line 



The search for stationary points of the error function in equation 22 (including the minimum point) leads to the 
two equations: 

(A;^) sin[26l] - 2C^y cos[26l] = (Al) 

{x) sin [61] - (y) cos[6l] + c = . (A2) 

The solution of equation | A 1 1 uses standard and well-known trigonometric procedures. However some care is required 
for the proper handling of the trigonometric functions. After the angle 9 is found, equation |A2| allows to determine c 
in a trivial way. 

Equat ion I A 1 1 al way s has two distinct solutions for 29, differing by ±7r, in the interval < 26* < 27r, as it is well-known 
from elementary trigonometry. In fact, letting X = cos[26'] and Y = sin[2^?] and using +1"^ = 1 the interpretation 
of the equation in terms of analytic geometry immediately leads to the conclusion. 

Therefore there are four distinct solutions for 0, in the interval Q < 29 < 27r, differing by ±7r/2. Two of them, 
differing by ±7r, minimize the error function and corresponds to the same straight line. The other two solutions, 
orthogonal to the first ones, maximize the error function. 

First note that the sign of the product C^y is the same as the sign of the product sin[2^?] cos [26*], so that this 
determines in which quadrant the angle 29 lies: 

l^C^y^Q <^ sin [261] cos [261] ^ . (A3) 

Second, the following special cases arise: 

Ay = <^ cos[26i] = 9^^+q'^ qeZ (A4) 

C^y^O ^ am[29] = ^ ^ ^ ^\ ^^'^ • ^^^^ 



Finally, if Ay ^ equation XT can be safely reduced to: 

Al/ ^ ^ tan[2^?] - ^ = ^ arctan[ ^] + q^- qeZ . (A6) 



Equation A6 for 9 provides two angles in the range — 7r/2 < 9 < +7r/2, differing by n/2, one corresponding to a 
minimum and the other to a maximum of the error function. 
Equation |A6| immediately implies the useful relations: 

cos2[20] = ^-f^ sm^[29] = . (A7) 

iAVf+ ACly {AVf + 4C^,y 



In order to determine which solution in equation A6 corresponds to the minimum/maximum it is easiest to re-start 
from equation |22[ and re- write it as follows, using a few trigonometric transformations: 

/c=l 

1 ^ 

^ nT. ((^fe - - iVk - {y))cos[9]f = (A9) 

fe=i 

= sin^ [9] + Vy cos^ [9] - 2C^y sin[6l] cos[6'] = (AlO) 
= \ ((Ke + Vy) - AF cos[26i] - 2C^y sin[26l]) = (All) 

1 / (Al/)^ + 4C2 \ 

= 2 (V^. + Vy) - cos[20] = (for AV ^ 0) (A12) 



\Uy.^ Vy) - sin[20] ^^^2cl^^"' ) = (for C,y ^ 0) (A13) 



^ {{Va, + Vy) ~ (AF)^ + ACly ) at the minimmn . (A14) 
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The above expressions clearly show that the minimum is obtained for those values of such that AFcos[26'] > 
and Cxy sm[20\ > 0. Therefore: 

AV 

Aycos[26i]>0 =^ cos [261] = — , (A15) 

2C 

C^y sm[20] > =^ sin[26'] = "'^ = . (A16) 



The above equations show that the angle 26 corresponding to the best-fit line lies in the first/fourth quadrant 
according to the positive/negative sign of C^y 

Clearly, the minimum value of the error function is zero for a perfect line fit: C^y = VxVy. 

The minimum of the error function can finally be expressed, in terms of a = Sign[C2;j,] = C^y/ \Cxy\ which makes 
sure that sin[26'] has the correct sign, as: 

a ( /sy 

arccos — 1^^^=^^^= \ + qtt qeZ choosing \9\ < tt/2 . (A17) 



2 



Alternatively the minimum/maximum can be determined as follows. Taking the second derivative with respect to 
of the equation |A10[ after using the stationary point conditions for both c and 9, one finds: 

f T^X^\0]\ , , , , sin[26'] /, o „N cos [261] . „x ^ ^ 

^ ( = ^ '"'[^^l + ^^^^ = ^ V^^^ + J = ^ V^^^ + ^^'y) ' ^^^^^ 

showing again that the minimum is given by those values of such that sin [261] and Cxy have the same sign. 



Appendix B: Limit of the OLS-y:x/OLS-x:y fit 



The limiting cases of the OLS-y:x/OLS-x:y fit, that is the case when one of the two variables has negligible errors, 
can be recovered. Consider a fixed set of N measured data points, {xk ± cj; yk ± (Jy} ,{k — 1, . . . , N), and let us 
study the limiting case when one of the standard errors tends to zero. 

The relation between the raw and re-scaled slopes is: 

tan[^] = tan[0] . (Bl) 

CTs: 



1. Limit for tlie slope/intercept 

Consider, for definiteness, the case OLS-y:x. 
Consider first the relations for the slope. 

In the limit — > 0, for equation [TOj one finds, the expression for the slope of the OLS-y:x fit as follows: 

Py ^ tan[05] ^^(Ay + aJl + AA Ay ee (B2) 



Consider now the relations for the intercepts: the equation in 10 for q is the same as the one for OLS-y:x fit, so 
there is need to investigate further its limit. 

Similarly, one may proceed for the limiting case (jy — > 0, to find the OLS-x:y fit limit. 
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2. Limit for the standard errors 



First re- write equation |29| in terms of the raw variables: 



N , 



X y xy 



vr 
2 



The identification of the angle 9 with the angles in equation [7] (see also equation 12) is obviously: 

69 

The relation between the error on the raw angle and the re-scaled angle is then: 

+ cr2 sin^fei 
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oy (l+ tan2 [9] \^^_cjI + cjI tan^ [9] 1 ^ a| cos^ [g] + gj sin ^ [9] 



<7x \ l+ tan^[9] J 1 + tan2[6'] (^xC^y 

From the above relation one can find for the raw angle: 



CTxCTv 



59. 



■cos'^[9] 59 
sm^^]59 



The above relation gives the error on the raw slope: 

c — > 5py = ^f + tan-^[6yj j uuy 

^ — ^ =^ 5p^={l+ tan2[a,] 



<7x 
Ox 
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N ^jVx 

1 crj 



f7. 



iV ^JVy 



(B4) 



(B5) 



(B6) 



(B7) 



(B8) 



exactly the ones for the OLS-y:x/OLS-x:y fit. 

Again, as the equations for the intercepts are the same as for the OLS-y:x/OLS-x:y fit, there is no need to investigate 
the limiting case. 



Appendix C: Derivation of the standard errors for slope and intercept 



The errors on the slope and intercept (see section III ) can be easily calculated as a by-product of the results obtained 
for the angle/signed-distance parametrization, starting from: 

p\9, c] = py = tan[6l] q\9, c] = qy = c7l + tan2[6'] = c/ cos[6'] cos[6'] > , (CI) 



without neglecting the covariance term in equation 31 Note that the sign of c is the same as the sign of qy 



{5vf ^{\-rV^)\59f 



{Sqf = (1 ( ^ + {59f [Z - cpf^ =^i+p^)^^+ i^59f {x)^ (1 



Cov[p, q] = - (1 (z - cp) {59f = - (x) i59y 



(C2) 
(C3) 
(C4) 



Equation 



C3 



has a simple interpretation by observing that (Z — cpY 



[l + p^) is the squared-distance between 



the centroid of the measured data points and the point where the best-fit line crosses the y axis. As the best-fit line 
always passes by the centroid, the term (l + p"^) (59)'^ (Z — cp)^ is the uncertainty on the intercept caused by the 
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uncertainty of the angle, {S9) . It is summed in quadrature to the term giving the error on the location of the 
centroid. 

It can be easily shown, by using the results of appendix [B| in particular equation |B6[ that the above formulas |C2[|C3| 
and |C4[ tend to the well-known results for the OLS fit. 
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